/*
 * HkfModelTest.cpp
 *
 *  Created on: 16 Aug 2011
 *      Author: Allan
 */

#include "HkfModelTest.h"

// C++ includes
#include <iostream>
#include <iomanip>
#include <cmath>

// Boost includes
#include <boost/foreach.hpp>

// GeoReact includes
#include "../../Applications/Equilibrator.h"
#include "../../Utils.h"

void TestHkfModel_WaterNaCl(double T, double P)
{
	Equilibrator equilibrator({"H2O", "NaCl(aq)"});
	
	cout << setw(15) << left << "mNaCl";
	cout << setw(15) << left << "ionic-strength";
	cout << setw(15) << left << "act-water";
	cout << setw(15) << left << "act-coeff-NaCl";
	cout << endl;

	BOOST_FOREACH(double mNaCl, Range(0.0, 6.0, 0.1))
	{
		equilibrator.ClearConditions();
		
		equilibrator.AddCondition("charge-balanced");
		equilibrator.AddCondition("concentration", "H2O", 55.508);
		equilibrator.AddCondition("total-concentration", "Na[+]", mNaCl);
		equilibrator.AddCondition("total-concentration", "Cl[-]", mNaCl);
		
		equilibrator.Solve(T, P);
		
		const double gNa = equilibrator.Get("activity-coeff", "Na[+]");
		const double gCl = equilibrator.Get("activity-coeff", "Cl[-]");
		
		const double gNaCl = std::sqrt(gNa * gCl);
		
		cout << setw(15) << left << mNaCl;
		cout << setw(15) << left << equilibrator.Get("ionic-strength");
		cout << setw(15) << left << equilibrator.Get("activity", "H2O");
		cout << setw(15) << left << gNaCl;
		cout << endl;
	}
}

void TestHkfModel_WaterCaCl2(double T, double P)
{
	Equilibrator equilibrator({"H2O", "CaCl2(aq)"});
	
	cout << setw(15) << left << "mCaCl2";
	cout << setw(15) << left << "ionic-strength";
	cout << setw(15) << left << "act-water";
	cout << setw(15) << left << "act-coeff-CaCl2";
	cout << endl;
	
	BOOST_FOREACH(double mCaCl2, Range(0.0, 3.0, 0.1))
	{
		equilibrator.ClearConditions();
		
		equilibrator.AddCondition("charge-balanced");
		equilibrator.AddCondition("concentration", "H2O", 55.508);
		equilibrator.AddCondition("total-concentration", "Ca[2+]", mCaCl2);
		equilibrator.AddCondition("total-concentration", "Cl[-]", 2.0*mCaCl2);
		
		equilibrator.Solve(T, P);
		
		const double gCa = equilibrator.Get("activity-coeff", "Ca[2+]");
		const double gCl = equilibrator.Get("activity-coeff", "Cl[-]");
		
		const double gCaCl2 = std::pow(gCa * gCl * gCl, 1.0/3);
		
		cout << setw(15) << left << mCaCl2;
		cout << setw(15) << left << equilibrator.Get("ionic-strength");
		cout << setw(15) << left << equilibrator.Get("activity", "H2O");
		cout << setw(15) << left << gCaCl2;
		cout << endl;
	}
}
